
/** include predefined life command line options */
#include <life/options.hpp>

/** include linear algebra backend */
#include <life/lifealg/backend.hpp>

/** include function space class */
#include <life/lifediscr/functionspace.hpp>

/** include helper function to define \f$P_0\f$ functions associated with regions  */
#include <life/lifediscr/region.hpp>

/** include integration methods */
#include <life/lifepoly/im.hpp>
/** include gmsh mesh importer */
#include <life/lifefilters/gmsh.hpp>

/** include exporter factory class */
#include <life/lifefilters/exporter.hpp>

/** include gmsh generator for tensorized domains */
//#include <life/lifefilters/gmshtensorizeddomain.hpp>

/** include  polynomialset header */
#include <life/lifepoly/polynomialset.hpp>

/** include  the header for the variational formulation language (vf) aka FEEL++ */
#include <life/lifevf/vf.hpp>

#include <life/lifefilters/geotool.hpp>
//#include <fstream> 

/** use Life namespace */
using namespace Life;
using namespace Life::vf;


inline
po::options_description
makeOptions()
{
  po::options_description code1doptions("Code1D options");
  code1doptions.add_options()
    ("hsize", po::value<double>()->default_value( 0.1 ), "mesh size")
    ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
    ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
    ("nu", po::value<double>()->default_value( 0.5 ), "nu")
    ("alphaX", po::value<double>()->default_value( 0.01 ), "alphaX")
    ("alphaY", po::value<double>()->default_value( -2.0 ), "alphaY")
    ("bccoeff", Life::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
    ("toler", Life::po::value<double>()->default_value( 1e-6 ), "error tol")
    ("xc", Life::po::value<double>()->default_value( 0.0 ), "param de recouvrement des maillages")
    ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
    ;
  return code1doptions.add( Life::life_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Life::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
    AboutData about( "Scharzz" ,
                     "Scharzz" ,
                     "0.2",
                     "Code 1D for -Delta u = f",
                     Life::AboutData::License_GPL,
                     "Copyright (c) 2010 Universite Joseph Fourier");

    about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
    return about;

}


template<uint Dim,uint N,uint OrderGeo>
class Localisation
  :
  public Application
{
  typedef Application super;
public:

  //! Polynomial order \f$P_2\f$
  //    static const uint16_type Order = 5;
  ///static const uint Order = N;

  typedef double value_type;

  //---------------------------------------------------------------------------------------//

  typedef Backend<value_type> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  typedef backend_type::sparse_matrix_type sparse_matrix_type;
  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
  typedef backend_type::vector_type vector_type;
  typedef backend_type::vector_ptrtype vector_ptrtype;

  //---------------------------------------------------------------------------------------//

  typedef Simplex<Dim,OrderGeo,Dim> convex_type;
  typedef Mesh<convex_type> mesh_type;
  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

  //---------------------------------------------------------------------------------------//

  typedef bases<Lagrange<N,Scalar/*,PointSetFekete*/> > basis_type;
  typedef FunctionSpace<mesh_type, basis_type> space_type;
  typedef boost::shared_ptr<space_type> space_ptrtype;
  typedef typename space_type::element_type element_type;

  //---------------------------------------------------------------------------------------//

  typedef Exporter<mesh_type> export_type;
  typedef boost::shared_ptr<export_type> export_ptrtype;

    /**
     * Constructor
     */
  Localisation( int argc, char** argv, AboutData const& ad,
		po::options_description const& od )
    :
    super( argc, argv, ad, od ),
    M_backend( backend_type::build( this->vm() ) ),
    M_backend2( backend_type::build( this->vm() ) ),
    //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
    meshSize( this->vm()["hsize"].template as<double>() ),
    meshSize1( this->vm()["hsize1"].template as<double>() ),
    meshSize2( this->vm()["hsize2"].template as<double>() ),
    exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) ),
    exporter1( Exporter<mesh_type>::New( this->vm(), "domain_1" ) ),
    exporter2( Exporter<mesh_type>::New( this->vm(), "domain_2" ) )
  {
  }

  mesh_ptrtype createMesh( std::string __name,std::string __str) ;

  mesh_ptrtype createMesh( double meshSize );
mesh_ptrtype createMesh1( double meshSize ,double x_couvr);
mesh_ptrtype createMesh2( double meshSize, double x_couvr,std::string leostr);
mesh_ptrtype createMesh3( double meshSize, double x_couvr,std::string leostr);
    
  // run the application
  void run();
  void runMultiDomain();
void runTest();
private:

    // solve the system \f$D u = F\f$
  void solve( sparse_matrix_ptrtype const& D, element_type& u, vector_ptrtype const& F, bool is_sym );
  //void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
  void problem(element_type& u , element_type& v , sparse_matrix_ptrtype& A , vector_ptrtype& B );
  void problem1(element_type& u , element_type& v , element_type& u_autre  );
  void problem2(element_type& u , element_type& v , element_type& u_autre );

    // export results to ensight format (enabled by  --export cmd line options)
  void exportResults( element_type& u );
  void exportResults( element_type& u1,element_type& u2,double time );

private:

  // linear algebra backend
  backend_ptrtype M_backend,M_backend2;
  sparse_matrix_ptrtype A1;
  vector_ptrtype B1;
  sparse_matrix_ptrtype A2;
  vector_ptrtype B2;
  // mesh characteristic size
  double meshSize,meshSize1,meshSize2;

  // exporter factory
  export_ptrtype exporter;
  export_ptrtype exporter1;
  export_ptrtype exporter2;
}; // Interpolate



template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str) 
{
  mesh_ptrtype mesh( new mesh_type );

  Gmsh gmsh;
  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( __name, __str  );

  ImporterGmsh<mesh_type> import( fname );
    
  float gmsh_version = this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );

  return mesh;

}


template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( double meshSize ) {
    mesh_ptrtype mesh( new mesh_type );

    std::ostringstream ostr;
    if (true)
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << meshSize << ";\n"
         << "Point(1) = {-1,-1, 0, h};\n"
         << "Point(2) = { 0,-1, 0, h};\n"
         << "Point(3) = { 1,-1, 0, h};\n"
         << "Point(4) = { 1, 1, 0, h};\n"
         << "Point(5) = { 0, 1, 0, h};\n"
         << "Point(6) = {-1, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Line(2) = {2, 3};\n"
         << "Line(3) = {3, 4};\n"
         << "Line(4) = {4, 5};\n"
         << "Line(5) = {5, 6};\n"
         << "Line(6) = {6, 1};\n"
         << "Line(7) = {2, 5};\n"
         << "Line Loop(1) = {1, 7, 5 ,6 };\n"
         << "Line Loop(2) = {2, 3, 4 , -7};\n"
         << "Plane Surface(2) = {1};\n"
         << "Plane Surface(3) = {2};\n"
         << "Physical Line(\"Bord\") = {1,2,3,4,5,6};\n"
         << "Physical Line(\"Interface\") = {7};\n"
         << "Physical Surface(\"Physique1\") = {2};\n"
         << "Physical Surface(\"Physique2\") = {3};\n";

    Gmsh gmsh;
    gmsh.setOrder(OrderGeo);
    std::string fname = gmsh.generate( "domain", ostr.str()  );

    ImporterGmsh<mesh_type> import( fname );
    
    float gmsh_version = this->vm()["gmsh"].template as<float>();
    if (gmsh_version==2.0) import.setVersion( "2.0" );
    else if (gmsh_version==2.1) import.setVersion( "2.1" );

    mesh->accept( import );

    return mesh;


}

template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh1( double meshSize,double x_couvr )
{
    mesh_ptrtype mesh( new mesh_type );
    
    double esperance=0.4;

    std::ostringstream ostr;
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << meshSize << ";\n"
         << "a="<<x_couvr<<";\n"
         << "a="<<esperance<<";\n"
         << "Point(1) = {-1,-1, 0, h};\n"
         << "Point(2) = { pt+a,-1, 0, h};\n"
         << "Point(3) = { pt+a, 1, 0, h};\n"
         << "Point(4) = {-1, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Line(2) = {2, 3};\n"
         << "Line(3) = {3, 4};\n"
         << "Line(4) = {4, 1};\n"
         << "Line Loop(1) = {1, 2, 3 ,4 };\n"
         << "Plane Surface(2) = {1};\n"
         << "Physical Line(\"Bord\") = {1,3,4};\n"
         << "Physical Line(\"Interface\") = {2};\n"
         << "Physical Surface(\"Physique1\") = {2};\n";


    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % Order
                            % this->vm()["hsize"].template as<double>()
                            );

    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain1", ostr.str()  );

    ImporterGmsh<mesh_type> import( fname );
    import.setVersion( "2.1" );
    mesh->accept( import );

    return mesh;

}
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//

template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh2( double __h, double x_couvr,std::string leostr )
{
    mesh_ptrtype mesh( new mesh_type );

    double esperance=0.4;

    /*    std::ostringstream ostr;
    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << __h << ";\n"
         << "a="<<-x_couvr<<";\n"
	 << "pt="<<esperance<<";\n"
         << "Point(1) = { pt-2*a,-1, 0, h};\n"
         << "Point(2) = { 1,-1, 0, h};\n"
         << "Point(3) = { 1, 1, 0, h};\n"
         << "Point(4) = { pt-2*a, 1, 0, h};\n"
         << "Line(1) = {1, 2};\n"
         << "Line(2) = {2, 3};\n"
         << "Line(3) = {3, 4};\n"
         << "Line(4) = {4, 1};\n"
         << "Line Loop(1) = {1, 2, 3 ,4 };\n"
         << "Plane Surface(2) = {1};\n"
         << "Physical Line(\"Bord\") = {1,2,3};\n"
         << "Physical Line(\"Interface\") = {4};\n"
         << "Physical Surface(\"Physique1\") = {2};\n";
    */
    std::ostringstream ostr;
    /*  ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=" << __h << ";\n"
         << "a="<<-x_couvr<<";\n"
	 << "pt="<<esperance<<";\n"
         << "Point(1) = { 0.7, 0.5, 0, h};\n"
         << "Point(2) = { 0.5, 0.5, 0, h};\n"
         << "Point(3) = { 0.3,0.5, 0, h};\n"
         << "Point(4) = { 0, 0, 0, h};\n"
         << "Point(5) = { 1, 0, 0, h};\n"
         << "Point(6) = { 1, 1, 0, h};\n"
         << "Point(7) = { 0, 1, 0, h};\n"
         << "Line(3) = {4, 5};\n"
         << "Line(4) = {5, 6};\n"
         << "Line(5) = {6, 7};\n"
         << "Line(6) = {7, 4};\n"
         << "Circle(1) = {1,2,3};\n"
         << "Circle(2) = {3,2,1};\n"
         << "Line Loop(2) = {3, 4, 5 ,6 };\n"
         << "Line Loop(1) = {1,2};\n"
         << "Plane Surface(1) = {2,1};\n"
         << "Physical Line(\"Bord\") = {3,4,5,6};\n"
         << "Physical Line(\"Interface\") = {1,2};\n"
         << "Physical Surface(\"Physique1\") = {1};\n";
*/

    ostr << "Mesh.MshFileVersion = " << 2 << ";\n"
         << "h=0.05;\n"
         << "Point(1) = { 0.7, 0.5, 0, h};\n"
         << "Point(2) = { 0.5, 0.5, 0, h};\n"
         << "Point(3) = { 0.3,0.5, 0, h};\n"
         << "Circle(1) = {1,2,3};\n"
         << "Circle(2) = {3,2,1};\n"
         << "Line Loop(1) = {1,2};\n"
         << "h=0.05;\n"
         << "Point(4) = { 0, 0, 0, h};\n"
         << "Point(5) = { 1, 0, 0, h};\n"
         << "Point(6) = { 1, 1, 0, h};\n"
         << "Point(7) = { 0, 1, 0, h};\n"
         << "Line(3) = {4, 5};\n"
         << "Line(4) = {5, 6};\n"
         << "Line(5) = {6, 7};\n"
         << "Line(6) = {7, 4};\n"
         << "Line Loop(2) = {3, 4, 5 ,6 };\n"
         << "Plane Surface(1) = {2,1};\n";
      // << "Physical Line(\"Interface\") = {1,2};\n"
      // << "Physical Line(\"Bord\") = {3,4,5,6};\n"
      //         << "Physical Surface(\"Omega1\") = {1};\n";

    std::cout<<ostr.str();

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % Order
                            % this->vm()["hsize"].template as<double>()
                            );


    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain2", leostr  );
    //std::string fname = gmsh.generate( "domain2", R.geoStr()  );
    //std::string fname = gmsh.generate( "domain2", ostr2.str()  );

    ImporterGmsh<mesh_type> import( fname );
    import.setVersion( "2.1" );
    mesh->accept( import );

    return mesh;

}

template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh3( double __h, double x_couvr,std::string leostr )
{
    mesh_ptrtype mesh( new mesh_type );

    double esperance=0.4;


    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % Order
                            % this->vm()["hsize"].template as<double>()
                            );

    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain2", leostr  );
    //std::string fname = gmsh.generate( "domain2", R.geoStr()  );
    //std::string fname = gmsh.generate( "domain2", ostr2.str()  );

    ImporterGmsh<mesh_type> import( fname );
    import.setVersion( "2.1" );
    mesh->accept( import );

    return mesh;

}


template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::exportResults( element_type& U )
{
    Log() << "exportResults starts\n";
    exporter->step(0)->setMesh( U.functionSpace()->mesh() );
    exporter->step(0)->add( "u", U );
    exporter->save();

} // Interpolate::export

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::exportResults( element_type& u1, element_type& u2, double time  )
{
    exporter1->step( time )->setMesh( u1.functionSpace()->mesh() );
    exporter1->step( time )->add( "u1", u1 );
    exporter1->save();

    exporter2->step( time )->setMesh( u2.functionSpace()->mesh() );
    exporter2->step( time )->add( "u2", u2 );
    exporter2->save();

} // Test_Schwarz::exportResults

template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::run()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

    using namespace Life::vf;  
 
    //-----------------------------------------------------------------------------------//

  int MyOrder=N;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % MyOrder
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  mesh_ptrtype mesh = createMesh( meshSize );

  space_ptrtype Xh = space_type::New( mesh );

  element_type u( Xh, "u" );
  element_type v( Xh, "v" );

  //-----------------------------------------------------------------------------------//

  AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

  vector_ptrtype R( M_backend->newVector( Xh ) );
  sparse_matrix_ptrtype J( M_backend->newMatrix(Xh,Xh) );

  this->problem(u,v,J,R);

  exportResults(u);

}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::runMultiDomain()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

    using namespace Life::vf;  
 
    //-----------------------------------------------------------------------------------//

    //int MyOrder=N;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % Order
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//


  GeoTool::Node x1(0,0);
  GeoTool::Node x2(1,1);
  GeoTool::Rectangle R(meshSize1,"RIGHT",x1,x2);
  R.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
  R.setMarker(_type="surface",_name="Omega1",_marker1=true);
  
  GeoTool::Node x_centre(0.5,0.5);
  GeoTool::Node x_bord(0.7,0.5);
  GeoTool::Circle C(meshSize2,"uncercle",x_centre,x_bord);
  C.setMarker(_type="line",_name="Interface",_marker1=true);
  //C.setMarker(_type="surface",_name="Omega2",_marker1=true);

  GeoTool::Node x2_centre(0.5,0.5);
  GeoTool::Node x2_bord(0.8,0.5);
  GeoTool::Circle C2(meshSize2,"uncercle",x2_centre,x2_bord);
  C2.setMarker(_type="line",_name="Interface",_marker1=true);
  C2.setMarker(_type="surface",_name="Omega2",_marker1=true);


  
  std::ostringstream ostr;
  ostr << (R-C).geoStr();// << "Physical Surface(\"Omega1\") = {1};\n";
  std::cout << ostr.str();
  mesh_ptrtype mesh1 = createMesh("mesh1",ostr.str());


  std::ostringstream ostr2;
  ostr2 << C2.geoStr();// << "Physical Surface(\"Omega2\") = {1};\n";
  //std::cout << ostr2.str();
  mesh_ptrtype mesh2 = createMesh("mesh2",ostr2.str());


  /*
    GeoTool::Node x3(-1,-1);
    GeoTool::Node x4(0,1);
  GeoTool::Rectangle R2(meshSize1,"kokok",x3,x4);
  R2.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R2.setMarker(_type="line",_name="Interface",_marker2=true);

  GeoTool::Node x5(0,-1);
  GeoTool::Node x6(1,1);
  GeoTool::Rectangle R3(meshSize2,"sgqysg",x5,x6);
  R3.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true);
  R3.setMarker(_type="line",_name="Interface",_marker4=true);

  std::ostringstream ostr3;
  ostr3 << R2.geoStr() << "Physical Surface(\"Omega1\") = {1};\n";
  //mesh_ptrtype mesh1 = createMesh("mesh1",ostr3.str());


  std::ostringstream ostr4;
  ostr4 << R3.geoStr() << "Physical Surface(\"Omega2\") = {1};\n";
  //mesh_ptrtype mesh2 = createMesh("mesh2",ostr4.str());
  */
  //-----------------------------------------------------------------------------------//


  //double xc = this->vm()["xc"].template as<double>();
  //mesh_ptrtype mesh2 = createMesh1( meshSize1,xc );
  //mesh_ptrtype mesh1 = createMesh2( meshSize2,xc,ostr.str() );


  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh2 = space_type::New( mesh2 );

  element_type u1( Xh1, "u1" );
  element_type v1( Xh1, "v1" );
  element_type u1b( Xh1, "u1b" );
  element_type u2( Xh2, "u2" );
  element_type v2( Xh2, "v2" );
  element_type u2b( Xh2, "u2b" );


  double __ppp =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
  double __ppp2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
  
  //-----------------------------------------------------------------------------------//

  A1=M_backend->newMatrix(Xh1,Xh1);
  B1=M_backend->newVector( Xh1 );

  A2= M_backend->newMatrix(Xh2,Xh2);
  B2= M_backend->newVector( Xh2 );

  //-----------------------------------------------------------------------------------//

  AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );
  
  double erreur=1;
  double err_tol=this->vm()["toler"].template as<double>();
  double err1,err2,norm1,norm2;
  int cpt=0;
    
  while (erreur>=err_tol) {
   
    u1b=u1;

    this->problem1(u1,v1,u2);
    norm1 =math::sqrt( integrate( elements(mesh1), _Q<2*Order>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
    err1=math::sqrt( integrate(elements(mesh1), _Q<2*Order>(), (idv(u1)-idv(u1b))*(idv(u1)-idv(u1b)) ).evaluate()(0,0) );
    if (norm1!=0)  err1=err1/norm1;

    u2b=u2;
    this->problem2(u2,v2,u1);
    norm2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
    err2=math::sqrt( integrate(elements(mesh2), _Q<2*Order>(), (idv(u2)-idv(u2b))*(idv(u2)-idv(u2b)) ).evaluate()(0,0) );
    if (norm2!=0)  err2=err2/norm2;

    erreur=err1+err2;
    std::cout<<"\n erreur = "<<erreur<<" ; err1 = "<<err1<<" ; err2 = "<<err2;

    //    exportResults( u1, u2, cpt);
    ++cpt;

    }
  

  //AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

  //  u1= project( Xh1, elements(mesh1),  SolAnalytique );

  //u2= project( Xh2, elements(mesh2),  SolAnalytique );

        exportResults( u1, u2, 0);


}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::runTest()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

    using namespace Life::vf;  
 
    //-----------------------------------------------------------------------------------//

    //int MyOrder=N;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % Order
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//
  double xc = this->vm()["xc"].template as<double>();
  mesh_ptrtype mesh1 = createMesh1( meshSize,xc );
  mesh_ptrtype mesh2 = createMesh2( meshSize,xc );

  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh2 = space_type::New( mesh2 );

  element_type u1( Xh1, "u1" );
  element_type v1( Xh1, "v1" );
  element_type u1b( Xh1, "u1b" );
  element_type u1ex( Xh1, "u1ex" );
  element_type u2( Xh2, "u2" );
  element_type v2( Xh2, "v2" );
  element_type u2b( Xh2, "u2b" );
  element_type u2ex( Xh2, "u2ex" );

  //-----------------------------------------------------------------------------------//

  AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );

  u1ex = project( Xh1, elements(mesh1),  SolAnalytique );
  u2ex = project( Xh2, elements(mesh2),  SolAnalytique );

  sparse_matrix_ptrtype A( M_backend->newMatrix(Xh1,Xh1) );
  vector_ptrtype B( M_backend->newVector( Xh1 ) );
  // A1->zero();
  //B1->zero();


  value_type pi = M_PI;
  value_type nu=0.5;

  AUTO( Alpha, vec(cst(0.01),cst(-2.)) );
  AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
  AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

  //  integrate( elements(mesh1), _Q<2*Order>(), (gradv(u1ex)*trans(gradv(u1ex) ).evaluate()(0,0) );

  AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);
  //------------------------------------------//
  form1( Xh1, B, _init=true ) =
    integrate( elements(mesh1), _Q<7>(),
	       f*id(v1) );

  form1( Xh1, B ) +=
    integrate( markedfaces(mesh1,mesh1->markerName("Interface")), _Q<2*Order>(),
	       nu*(gradv(u1ex/*_autre*/)*N())*id(v1) );

  B->close();
  std::cout<<"\nnomrL2 B =" << B->l2Norm();
  //------------------------------------------//
  B->zero();
  form1( Xh1, B, _init=true ) =
    integrate( elements(mesh1), _Q<7>(),
	       f*id(v1) );

  form1( Xh1, B ) +=
    integrate( markedfaces(mesh1,mesh1->markerName("Interface")), _Q<2*Order>(),
	       nu*(gradv(u2ex)*N())*id(v1) );

  B->close();
  std::cout<<"\nnomrL2 B =" << B->l2Norm();

  /*double erreur=1;
  double err_tol=1e-6;
  double err1,err2;
  int cpt=0;
  std::cout<<"\n--1--\n";
    while (erreur>=err_tol) {
    
    u1b=u1;
  std::cout<<"\n--2--\n";  
    this->problem1(u1,v1,u2);
    err1=math::sqrt( integrate(elements(mesh1), _Q<2*Order>(), (idv(u1)-idv(u1b))*(idv(u1)-idv(u1b)) ).evaluate()(0,0) );
  std::cout<<"\n--3--\n";
    u2b=u2;
    this->problem2(u2,v2,u1);
    err2=math::sqrt( integrate(elements(mesh2), _Q<2*Order>(), (idv(u2)-idv(u2b))*(idv(u2)-idv(u2b)) ).evaluate()(0,0) );
  std::cout<<"\n--4--\n";
    erreur=err1+err2;
    std::cout<<"\n erreur = "<<erreur<<" ; err1 = "<<err1<<" ; err2 = "<<err2;

    exportResults( u1, u2, cpt);
    ++cpt;

    }*/

  //exportResults(u);

}

  
template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::problem(element_type& u , element_type& v , sparse_matrix_ptrtype& A , vector_ptrtype& B ) {

    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    space_ptrtype Xh=u.functionSpace();
    mesh_ptrtype mesh=Xh->mesh();

    //----------------------------------------------------------------------------------------------------------//

    value_type pi = M_PI;
    value_type nu=0.5;
    double penalbc = this->vm()["bccoeff"].template as<value_type>();
    AUTO( Alpha, vec(cst(0.01),cst(-2.)) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

    AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);

    //----------------------------------------------------------------------------------------------------------//
    
    form2( Xh, Xh, A, _init=true ) =
        integrate( elements(mesh), _Q<2*2>(),
                   nu*gradt(u)*trans(grad(v)) );

    form2( Xh, Xh, A) +=
        integrate( elements(mesh), _Q<2*2>(),
                   (gradt(u)*Alpha)*id(v) );

    form2( Xh, Xh, A) +=
        integrate( markedfaces(mesh,mesh->markerName("Bord")), _Q<2*2>(),
		   - penalbc*nu*(gradt(u)*N())*id(v)
		   + penalbc*idt(u)*id(v)/hFace() );


    //----------------------------------------------------------------------------------------------------------//

   form1( Xh, B, _init=true ) =
        integrate( elements(mesh), _Q<2*2>(),
                   f*id(v) );

    //----------------------------------------------------------------------------------------------------------//


    A->close();

#if 0
    form2( Xh, Xh, A ) +=
      //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
      on( markedfaces(mesh, mesh->markerName("Bord")) , u, B, cst(0.) ) ;
#endif

    this->solve(A,u,B,false);

}


//--//////--//////--//////--//////-------------------------------------------------------------------------------------------------//
//-//--//--//--//--//--//--//--//---------------------------------------------------------------------------------------//
////--//////--//////--//////--////------------------------------------------------------------------------------------------------//

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::problem1(element_type& u , element_type& v , element_type& u_autre ) {

    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    space_ptrtype Xh=u.functionSpace();
    mesh_ptrtype mesh=Xh->mesh();

    sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
    vector_ptrtype B( M_backend->newVector( Xh ) );
  //    A1->zero();
  // B1->zero();

    //----------------------------------------------------------------------------------------------------------//

    value_type pi = M_PI;
    value_type nu= this->vm()["nu"].template as<value_type>();
    value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
    value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

    AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

    AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);




    //----------------------------------------------------------------------------------------------------------//
    
    form2( Xh, Xh, A, _init=true ) =
        integrate( elements(mesh), _Q<2*Order>(),
                   nu*gradt(u)*trans(grad(v)) );

    form2( Xh, Xh, A) +=
        integrate( elements(mesh), _Q<2*Order>(),
                   (gradt(u)*Alpha)*id(v) );

    //----------------------------------------------------------------------------------------------------------//

    form1( Xh, B, _init=true ) =
        integrate( elements(mesh), _Q<7>(),
                   f*id(v) );


    form1( Xh, B ) +=
        integrate( markedfaces(mesh,mesh->markerName("Interface")), _Q<2*Order>(),
		   nu*(gradv(u_autre)*N())*id(v) );

    //----------------------------------------------------------------------------------------------------------//

      A->close();
        form2( Xh, Xh, A ) +=
      //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
      on( markedfaces(mesh, mesh->markerName("Bord")) , u, B, cst(0.) );
  

  this->solve(A,u,B,false);


}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::problem2(element_type& u , element_type& v , element_type& u_autre ) {

    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//

    space_ptrtype Xh=u.functionSpace();
    mesh_ptrtype mesh=Xh->mesh();

    sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );
    vector_ptrtype B( M_backend->newVector( Xh ) );
    //A->zero();
    //B2->zero();
    //----------------------------------------------------------------------------------------------------------//

    value_type pi = M_PI;
    value_type nu= this->vm()["nu"].template as<value_type>();
    value_type alphaX =  this->vm()["alphaX"].template as<value_type>();
    value_type alphaY =  this->vm()["alphaY"].template as<value_type>();

    AUTO( Alpha, vec(cst(alphaX),cst(alphaY)) );
    AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    AUTO( lap_Uex, -(pi*pi/2.)*cos(pi*Px()/2.)*cos(pi*Py()/2.) );

    AUTO( f, trans(Alpha)*grad_Uex -nu*lap_Uex);

    //----------------------------------------------------------------------------------------------------------//

    form2( Xh, Xh, A, _init=true ) =
        integrate( elements(mesh), _Q<2*Order>(),
                   nu*gradt(u)*trans(grad(v)) );

    form2( Xh, Xh, A) +=
        integrate( elements(mesh), _Q<2*Order>(),
                   (gradt(u)*Alpha)*id(v) );
    
    /*
    //weakbc for interface
    double penalbc = this->vm()["bccoeff"].template as<double>();
    form2( Xh, Xh, A) +=
        integrate( markedfaces(mesh, mesh->markerName("Interface")), _Q<2*Order>(),
                   - nu*(gradt(u)*N())*id(v)
		   + penalbc*idt(u)*id(v)/hFace() ); 
    */
    //----------------------------------------------------------------------------------------------------------//

   form1( Xh, B, _init=true ) =
        integrate( elements(mesh), _Q<7>(),
                   f*id(v) );
   /*
   //weakbc for interface
   form1( Xh, B ) +=
        integrate( markedfaces(mesh, mesh->markerName("Interface")), _Q<2*Order>(),
		   penalbc*idv(u_autre)*id(v)/hFace() );
   */ 
 //----------------------------------------------------------------------------------------------------------//
   
#if 0
    A->close();
    form2( Xh, Xh, A ) +=
      //on( boundaryfaces(mesh) , u, B, cst(0.) ) ;
      on( markedfaces(mesh, mesh->markerName("Bord")) , u, B, cst(0.) );
#endif

    A->close();
    form2( Xh, Xh, A ) +=
      on( markedfaces(mesh, mesh->markerName("Interface")) , u, B, idv(u_autre) );

    this->solve(A,u,B,false);
    //    vector_ptrtype U( M_backend/*2*/->newVector( u.functionSpace() ) );
    //    M_backend2->solve( A2, A2, U, B2, false );
    //u = *U;



}



template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::solve( sparse_matrix_ptrtype const& D,
                                   element_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{

  backend_ptrtype M_backend2(backend_type::build( this->vm() ));
  
  vector_ptrtype U( M_backend2->newVector( u.functionSpace() ) );
  M_backend2->solve( D, D, U, F, false );
  u = *U;
} // Test_Schwarz::solve
 







/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
  //dim, N, OrderGeo
  Localisation<2,2,1> moncode( argc, argv, makeAbout(), makeOptions() );

  //  moncode.run();
  moncode.runMultiDomain();
  //  moncode.runTest();
}



